Lake Elsinore image processing
library(knitr)
opts_chunk$set(warning = FALSE, message = FALSE, fig.path = 'figs/', dev.args = list(bg = 'transparent'), eval = T)
library(raster)
library(sf)
library(tidyverse)
library(mapview)
library(proj4shortcut)
library(plotly)
prj_geo <- geo_wgs84
prj_pro <- utm_wgs84('11s')August 21st
Map of in situ samples with chlorophyll concentration.
# get coverage of in situ points by date
insit <- read_csv('raw/CHL_2017_08_21.csv', col_names = F) %>%
rename(
site = X1,
lat = X2,
lon = X3,
chl = X4
) %>%
mutate(lat = ifelse(site == 34, 33.66954, lat))
coordinates(insit) <- ~ lon + lat
proj4string(insit) <- prj_geo
mapview(insit, zcol = 'chl', legend = T) Overlap of in situ samples (red numbers) with the spatial extent of each image. Note that the spatial extent is rectangular and not all pixels in the extent have values.
# transform insit to projected
insit <- spTransform(insit, CRS(prj_pro))
# raster tiffs
imgs_pth <- list.files('img/GeoTiff_08_21/', pattern = '^Raster*', full.names = T)
# get raster extents
out <- list()
cmbs <- for(i in 1:length(imgs_pth)){
# raster aggregate, dissolve to polygon
rst_chk <- stack(imgs_pth[[i]]) %>%
.[[1]] %>%
aggregate(fact = 10, fun = mean)
values(rst_chk) <- ifelse(is.na(values(rst_chk)), NA, 1)
ply <- rst_chk %>%
rasterToPolygons(dissolve = T)
out[[i]] <- ply
}
# fortify polygons for plot
allply <- out %>%
map(., function(x){
dat <- fortify(x) %>%
mutate(fl = gsub('\\.[0-9]$', '', names(x@data)))
return(dat)
}) %>%
enframe %>%
unnest
# plot
p <- ggplot() +
geom_polygon(data = allply, aes(x = long, y = lat, group = fl), alpha = 0.5) +
geom_text(data = data.frame(insit), aes(x = lon, y = lat, label = site)) +
coord_equal()
ggplotly(p)Extracted pixel values for each of four bands in each image where sample sites occurred within an image.
# extract values from raster
out <- list()
for(i in 1:length(imgs_pth)){
# cat(i, 'of', length(imgs_pth), '\n')
# extract raster cells by points in insit
rst <- stack(imgs_pth[[i]])
crs(rst) <- prj_pro
rst_chk <- rst %>%
raster::extract(insit, buffer = 0) %>%
enframe %>%
bind_cols(insit@data, .) %>%
select(-name) %>%
filter(map(.$value, function(x) !anyNA(x)) %>% unlist) %>%
mutate(value = map(value, function(x){
sumpt <- x %>%
data.frame(din = .) %>%
rownames_to_column('img_bnd') %>%
filter(!grepl('\\.5$', img_bnd))
return(sumpt)
})) %>%
unnest
# append to output
out[[i]] <- rst_chk
}
# combine all extracted samples
# get ndvi, gndvi, vigreen
# current not in nm
exts <- out %>%
do.call('rbind', .) %>%
remove_rownames() %>%
separate(img_bnd, c('img', 'bnd'), sep = '\\.') %>%
select(-img) %>%
mutate(bnd = factor(bnd, levels = c('1', '2', '3', '4'), labels = c('grn', 'red', 'edg', 'nir'))) %>%
spread(bnd, din) %>%
mutate(
ndvi = (nir - red) / (nir + red), # normalized diff veg index
gndvi = (nir - grn) / (nir + grn), # green normalized diff veg index
vigrn = (grn - red) / (grn + red) # normalized diff of green and red bands
) %>%
gather('var', 'val', -site, -chl)
exts## # A tibble: 91 x 4
## site chl var val
## <int> <dbl> <chr> <dbl>
## 1 1 146 grn 0.267
## 2 2 143 grn 0.271
## 3 16 52.6 grn 0.269
## 4 17 209 grn 0.308
## 5 31 189 grn 0.330
## 6 32 157 grn 0.315
## 7 33 185 grn 0.337
## 8 34 170 grn 0.304
## 9 35 186 grn 0.302
## 10 36 285 grn 0.647
## # ... with 81 more rowsScatterplots of band values (din) with measured chlorophyll.
ggplot(exts, aes(x = val, y = chl)) +
geom_point() +
facet_wrap(~var, ncol = 2, scales = 'free_x') +
stat_smooth(method = 'lm') +
theme_bw()
Working with model selection:
library(caret)
library(MuMIn)
# https://pix4d.com/sequoia-faq/#3
tomod <- exts %>%
mutate(bnd = paste0('b', bnd)) %>%
spread(bnd, din)
glb <- lm(chl ~ (b1 + b2 + b2 + b3 + b4)^3 + I(b1^2) + I(b2^2) + I(b3^2) + I(b4^2), tomod,
na.action = 'na.pass')
tmp <- dredge(glb, evaluate = F)
# createFolds(1:nrow(tomod), k = 5)September 6th
Map of in situ samples with chlorophyll concentration.
# get coverage of in situ points by date
insit <- read_csv('raw/CHL_2017_09_06.csv', col_names = F) %>%
rename(
site = X1,
lat = X2,
lon = X3,
chl = X4
) %>%
mutate(lon = ifelse(site == 39, -117.3597, lon))
coordinates(insit) <- ~ lon + lat
proj4string(insit) <- prj_geo
mapview(insit, zcol = 'chl', legend = T) Overlap of in situ samples (red numbers) with the spatial extent of each image. Note that the spatial extent is rectangular and not all pixels in the extent have values.
# transform insit to projected
insit <- spTransform(insit, CRS(prj_pro))
# raster tiffs
imgs_pth <- list.files('img/GeoTiff_09_06/', pattern = '*\\.tif$', full.names = T)
# get raster extents
out <- list()
cmbs <- for(i in 1:length(imgs_pth)){
# raster aggregate, dissolve to polygon
rst_chk <- stack(imgs_pth[[i]]) %>%
.[[1]] %>%
aggregate(fact = 10, fun = mean)
values(rst_chk) <- ifelse(is.na(values(rst_chk)), NA, 1)
ply <- rst_chk %>%
rasterToPolygons(dissolve = T)
out[[i]] <- ply
}
# fortify polygons for plot
allply <- out %>%
map(., function(x){
dat <- fortify(x) %>%
mutate(fl = gsub('\\.[0-9]$', '', names(x@data)))
return(dat)
}) %>%
enframe %>%
unnest
# plot
p <- ggplot() +
geom_polygon(data = allply, aes(x = long, y = lat, group = fl), alpha = 0.5) +
geom_text(data = data.frame(insit), aes(x = lon, y = lat, label = site)) +
coord_equal()
ggplotly(p)Extracted pixel values for each of four bands in each image where sample sites occurred within an image.
# extract values from raster
out <- list()
for(i in 1:length(imgs_pth)){
# cat(i, 'of', length(imgs_pth), '\n')
rst <- stack(imgs_pth[[i]])
crs(rst) <- prj_pro
rst_chk <- rst %>%
raster::extract(insit) %>%
data.frame(insit@data, .) %>%
gather('img_bnd', 'din', -site, -chl) %>%
na.omit
# out <- c(out, list(tmp))
out[[i]] <- rst_chk
}
exts <- out %>%
do.call('rbind', .) %>%
remove_rownames() %>%
separate(img_bnd, c('img', 'bnd'), sep = '\\.') %>%
group_by(site, chl, bnd) %>%
summarise(din = mean(din)) %>%
ungroup %>%
mutate(bnd = factor(bnd, levels = c('1', '2', '3', '4'), labels = c('grn', 'red', 'edg', 'nir'))) %>%
spread(bnd, din) %>%
mutate(
ndvi = (nir - red) / (nir + red), # normalized diff veg index
gndvi = (nir - grn) / (nir + grn), # green normalized diff veg index
vigrn = (grn - red) / (grn + red) # normalized diff of green and red bands
) %>%
gather('var', 'val', -site, -chl)
exts## # A tibble: 70 x 4
## site chl var val
## <int> <dbl> <chr> <dbl>
## 1 1 62.6 grn 0.319
## 2 2 487 grn 0.322
## 3 3 185 grn 0.315
## 4 16 380 grn 0.326
## 5 17 295 grn 0.298
## 6 18 306 grn 0.305
## 7 19 293 grn 0.313
## 8 31 254 grn 0.384
## 9 32 215 grn 0.338
## 10 33 152 grn 0.497
## # ... with 60 more rowsScatterplots of band values (din) with measured chlorophyll.
ggplot(exts, aes(x = val, y = chl)) +
geom_point() +
facet_wrap(~var, ncol = 2, scales = 'free_x') +
stat_smooth(method = 'lm') +
theme_bw()